home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_200 / 228_01 / mathmax.c < prev    next >
Encoding:
Text File  |  1987-07-29  |  10.9 KB  |  505 lines
  1. /*
  2. HEADER:         CUGXXX;
  3. TITLE:          Mathematical subroutines;
  4. DATE:           3-20-86;
  5. DESCRIPTION:    Mathematical subroutines;
  6. KEYWORDS:       Mathematics, Math functions;  
  7. FILENAME:       MATHMAX.C;
  8. WARNINGS:       None;
  9. AUTHORS:        Max R. Dürsteler;
  10. COMPILER:       Lattice C, Microsoft C;
  11. REFERENCES:     US-DISK 1307;
  12. ENDREF
  13. */
  14.  
  15. /*****************************************************************************
  16.  *                                         *
  17.  * MATH.C  Mathematical subroutines for use with Lattice or Microsoft-         *
  18.  *       compiler.                                 *
  19.  *       By Max R. Dürsteler, 12405 Village Sq. Terrace, Rockville MD      *
  20.  *                                         *
  21.  *****************************************************************************
  22.  */
  23.  
  24. /* SQRT.C
  25.  * Calculates square root
  26.  * Precision: 11.05 E-10
  27.  * Range: 0 to e308
  28.  * Author: Max R. Dürsteler
  29.  */
  30. double sqrt(x)
  31. double x;
  32. {
  33.     typedef union {           /* Makes use of internal data  */
  34.         double d;          /* representation */
  35.         unsigned u[4];
  36.     } DBL;
  37.     double y, re, p, q;
  38.     unsigned ex;        /* exponens*/
  39.     DBL *xp;
  40.  
  41.     xp = (DBL *)&x;
  42.     ex = xp->u[3] & ~0100017; /* save exponens */
  43.     re = ex & 020 ? 1.4142135623730950488 : 1.0;
  44.     ex = ex - 037740 >> 1;
  45.     ex &= ~0100017;
  46.     xp->u[3] &= ~0177760;      /* erase exponens and sign*/
  47.     xp->u[3] |= 037740;      /* arrange for mantissa in range 0.5 to 1.0 */
  48.     if (xp->d < .7071067812) {  /* multiply by sqrt(2) if mantissa < 0.7 */
  49.         xp->d *= 1.4142135623730950488;
  50.         if (re > 1.0) re = 1.189207115002721;
  51.         else re = 1.0 / 1.18920711500271;
  52.     }
  53.     p =           .54525387389085 +   /* Polynomial approximation */
  54.         xp->d * (13.65944682358639 +   /* from: COMPUTER APPROXIMATIONS*/
  55.         xp->d * (27.090122712655   +   /* Hart, J.F. et al. 1968 */
  56.         xp->d *   6.43107724139784));
  57.     q =          4.17003771413707 +
  58.         xp->d * (24.84175210765124 +
  59.         xp->d * (17.71411083016702 +
  60.         xp->d ));
  61.     xp->d = p / q;
  62.     xp->u[3] = xp->u[3] + ex & ~0100000;
  63.     xp->d *= re;
  64.     return(xp->d);
  65. }
  66.  
  67.  
  68. /* EXP2.C
  69.  * Calculates exponens to base 2 using polynomial approximations
  70.  * Precision: 10E-10
  71.  * Range:
  72.  * Author: Max R. Dürsteler
  73.  */
  74. double exp2(x)
  75. double x;
  76. {
  77.     typedef union {
  78.         double d;
  79.         unsigned u[4];
  80.     } DBL;
  81.     double y, x2, p, q, re;
  82.     int ix;
  83.     DBL *xp, *yp;
  84.  
  85.     xp = (DBL *)&x;
  86.     y = 0.0;
  87.     yp = (DBL *)&y;
  88.     if(xp->d > 1023.0 || xp->d < -1023.0) return (1E307);
  89.     ix = (int) xp->d;
  90.     yp->u[3] += ix + 1023 << 4;
  91.     yp->u[3] &= ~0100017;
  92.     if ((xp->d -= (double) ix) == 0.0) return(yp->d);
  93.     if (xp->d < 0.0) {
  94.           yp->u[3] -= 1 << 4;
  95.           yp->u[3] &= ~0100017;
  96.           xp->d++;
  97.     }
  98.     if (xp->d >= 0.5) {  /* adjust to range 0-0.5 */
  99.         xp->d -= 0.5;
  100.         re = 1.41421356237309504880;
  101.     }
  102.     else re = 1.0;
  103.     x2 = xp->d * xp->d;
  104.     p = xp->d * (7.2152891511493 +
  105.            x2 *  0.0576900723731);
  106.     q =        20.8189237930062 +
  107.            x2 ;
  108.     xp->d = (q + p) / (q - p);
  109.     yp->d *= re * xp->d;
  110.     return(yp->d);
  111. }
  112.  
  113.  
  114. /* POW.C
  115.  * Calculates y^x
  116.  * Precision:
  117.  * Range: o to big for y, +/-1023 for x
  118.  * Author: Max R. Dürsteler, 10/2/83
  119.  */
  120.  
  121. double pow(y,x)
  122. double y, x;
  123. {
  124.     typedef union {
  125.         double d;
  126.         unsigned u[4];
  127.     } DBL;
  128.     double z, w, p, p2, q, re;
  129.     unsigned ex;        /* exponens*/
  130.     int iz;
  131.     DBL *yp, *zp, *wp;
  132.  
  133.     yp = (DBL *)&y;
  134.     if (yp->d <= 0.0) y = -y;
  135.     z = 0.0;
  136.     zp = (DBL *)&z;
  137.     zp->u[3] = yp->u[3] & ~0100017; /* save exponens */
  138.     iz = (zp->u[3] >> 4)-1023;
  139.     if ((yp->d - zp->d) == 0.0)
  140.         z = (double)iz;
  141.     else {
  142.         yp->u[3] -= ++iz << 4; /* arrange for range 0.5 to 0.99999999999 */
  143.         yp->d *= 1.4142135623730950488;  /* shift for 1/sqrt(2) to sqrt(2) */
  144.         p = (yp->d - 1.0) / (yp->d + 1.0);
  145.         p2 = p * p;
  146.         z = p  * (2.000000000046727 +    /* Polynomial approximation */
  147.             p2 * (0.666666635059382 +    /* from: COMPUTER APPROXIMATIONS*/
  148.             p2 * (0.4000059794795   +    /* Hart, J.F. et al. 1968 */
  149.             p2 * (0.28525381498     +
  150.             p2 *  0.2376245609 ))));
  151.         z = z * 1.442695040888634 + (double)iz    - 0.5;
  152.     }
  153.     z *= x;
  154.     w = 0.0;
  155.     wp = (DBL *)&w;
  156.     if(zp->d > 1023.0 || zp->d < -1023.0) return (1E307);
  157.     iz = (int) zp->d;
  158.     wp->u[3] += iz + 1023 << 4;
  159.     wp->u[3] &= ~0100017;
  160.     if ((zp->d -= (double) iz) == 0.0) return(wp->d);
  161.     while (zp->d < 0.0) {
  162.           wp->u[3] -= 1 << 4;
  163.           wp->u[3] &= ~0100017;
  164.           zp->d++;
  165.     }
  166.     if (zp->d >= 0.5) {  /* adjust to range 0-0.5 */
  167.         zp->d -= 0.5;
  168.         re = 1.41421356237309504880;
  169.     }
  170.     else re = 1.0;
  171.     p2 = zp->d * zp->d;
  172.     p = zp->d * (7.2152891511493 +
  173.            p2 *  0.0576900723731);
  174.     q =        20.8189237930062 +
  175.            p2 ;
  176.     zp->d = re * wp->d * (q + p) / (q - p);
  177.     return(zp->d);
  178. }
  179.  
  180.  
  181. /* FMOD.C
  182.  * Returns the number f such that x = i*y + f. i is an integer, and
  183.  * 0 <= f < y.
  184.  */
  185. double fmod(x, y)
  186. double x, y;
  187. {
  188.     double zint, z;
  189.     int i;
  190.  
  191.     z = x / y;
  192.     zint = 0.0;
  193.     while (z > 32768.0) {
  194.         zint += 32768.0;
  195.         z -= 32768;
  196.     }
  197.     while (z < -32768.0) {
  198.         zint -= 32768.0;
  199.         z += 32768.0;
  200.     }
  201.     i = (int) z;
  202.     zint += (double) i;
  203.     return( x - zint * y);
  204. }
  205.  
  206. /*ASIN.C
  207.  *Calculates arcsin(x)
  208.  *Range: 0 <= x <= 1
  209.  *Precision: +/- .000,000,02
  210.  *Header: math.h
  211.  *Author: Max R. Dürsteler
  212.  */
  213. extern double sqrt();
  214.  
  215. double asin(x)
  216. double x;
  217.  
  218. {
  219.     double y;
  220.     int sign;
  221.  
  222.     if (x > 1.0 || x < -1.0) exit(1);
  223.     sign = 0;
  224.     if (x < 0) {
  225.          sign = 1;
  226.          x = -x;
  227.     }
  228.     y = ((((((-.0012624911    * x
  229.          + .0066700901) * x
  230.          - .0170881256) * x
  231.          + .0308918810) * x
  232.          - .0501743046) * x
  233.          + .0889789874) * x
  234.          - .2145988016) * x
  235.          +1.5707963050;
  236.     y = 1.57079632679 - sqrt(1.0 - x) * y;
  237.     if (sign) y = -y;
  238.     return(y);
  239. }
  240.  
  241. /* SIN.C
  242.  * Calculates sin(x), angle x must be in rad.
  243.  * Range: -pi/2 <= x <= pi/2
  244.  * Precision: +/- .000,000,005
  245.  * Header: math.h
  246.  * Author: Max R. Dürsteler
  247.  */
  248.  
  249. double sin(x)
  250. double x;
  251. {
  252.     double xi, y, q, q2;
  253.     int sign;
  254.  
  255.     xi = x; sign = 1;
  256.     while (xi < -1.57079632679489661923) xi += 6.28318530717958647692;
  257.     while (xi > 4.71238898038468985769) xi -= 6.28318530717958647692;
  258.     if (xi > 1.57079632679489661923) {
  259.         xi -= 3.141592653589793238462643;
  260.         sign = -1;
  261.     }
  262.     q = xi / 1.57079632679; q2 = q * q;
  263.     y = ((((.00015148419  * q2
  264.           - .00467376557) * q2
  265.           + .07968967928) * q2
  266.           - .64596371106) * q2
  267.           +1.57079631847) * q;
  268.     return(sign < 0? -y : y);
  269. }
  270.  
  271.  
  272. /* LOG.C
  273.  * Calculates natural logarithmus
  274.  * Precision: 11.56 E-10
  275.  * Range: 0 to e308
  276.  * Author: Max R. Dürsteler
  277.  */
  278. double log(x)
  279. double x;
  280. {
  281.     typedef union {
  282.         double d;
  283.         unsigned u[4];
  284.     } DBL;
  285.     double y, z, z2, p;
  286.     unsigned ex;        /* exponens*/
  287.     int ix;
  288.     DBL *xp, *yp;
  289.  
  290.     xp = (DBL *)&x;
  291.     if (xp->d <= 0.0) return(y);
  292.     y = 0.0;
  293.     yp = (DBL *)&y;
  294.     yp->u[3] = xp->u[3] & ~0100017; /* save exponens */
  295.     ix = (yp->u[3] >> 4)-1023;
  296.     if ((xp->d - yp->d) == 0.0) return( .693147180559945 * (double)ix);
  297.     xp->u[3] -= ++ix << 4; /* arrange for range 0.5 to 0.99999999999 */
  298.     xp->d *= 1.4142135623730950488;  /* shift for 1/sqrt(2) to sqrt(2) */
  299.     z = (xp->d - 1.0) / (xp->d + 1.0);
  300.     z2 = z * z;
  301.     y = z  * (2.000000000046727 +    /* Polynomial approximation */
  302.         z2 * (0.666666635059382 +    /* from: COMPUTER APPROXIMATIONS*/
  303.         z2 * (0.4000059794795   +    /* Hart, J.F. et al. 1968 */
  304.         z2 * (0.28525381498     +
  305.         z2 *  0.2376245609 ))));
  306.     y = y + .693147180559945 * ((double)ix    - 0.5);
  307.     return(yp->d);
  308. }
  309.  
  310. /* ATAN.C
  311.  * Calculates arctan(x)
  312.  * Range: -infinite <= x <= infinite (Output -pi/2 to +pi/2)
  313.  * Precision: +/- .000,000,04
  314.  * Header: math.h
  315.  * Author: Max R. Dürsteler 9/15/83
  316.  */
  317.  
  318. double atan(x)
  319. double x;
  320. {
  321.     double xi, q, q2, y;
  322.     int sign;
  323.  
  324.     xi = (x < 0. ? -x : x);
  325.     q = (xi - 1.0) / (xi + 1.0); q2 = q * q;
  326.     y = ((((((( - .0040540580 * q2
  327.              + .0218612286) * q2
  328.              - .0559098861) * q2
  329.              + .0964200441) * q2
  330.              - .1390853351) * q2
  331.              + .1994653599) * q2
  332.              - .3332985605) * q2
  333.              + .9999993329) * q + 0.785398163397;
  334.     return(x < 0. ? -y: y);
  335. }
  336.  
  337.  
  338. /* FLOOR.C
  339.  * Returns largest integer not greater than x
  340.  * Author: Max R. Dürsteler 9/26/83
  341.  */
  342. double floor(x)
  343. double x;
  344. {
  345.     double y;
  346.     int ix;
  347.  
  348.     y = 0.0;
  349.     while (x >= 32768.0) {
  350.         y += 32768.0;
  351.         x -= 32768.0;
  352.     }
  353.     while (x <= -32768.0) {
  354.         y -= 32768.0;
  355.         x += 32768.0;
  356.     }
  357.     if (x > 0.0) ix = (int) x;
  358.     else ix = (int)(x - 0.9999999999);
  359.     return( y + (double) ix);
  360. }
  361.  
  362.  
  363. /* CEIL.C
  364.  * Returns smallest integer not less than x
  365.  * Author: Max R. Dürsteler 9/26/83
  366.  */
  367. double ceil(x)
  368. double x;
  369. {
  370.     double y;
  371.     int ix;
  372.  
  373.     y = 0.0;
  374.     while (x >= 32768.0) {
  375.         y += 32768.0;
  376.         x -= 32768.0;
  377.     }
  378.     while (x <= -32768.0) {
  379.         y -= 32768.0;
  380.         x += 32768.0;
  381.     }
  382.     if (x > 0.0) ix = (int) (x + 0.999999999999999);
  383.     else ix = (int) x;
  384.     return( y + (double) ix);
  385. }
  386.  
  387. /* EXP.C
  388.  * Calculates exponens of x to base e
  389.  * Range: +/- exp(88)
  390.  * Precision: +/- .000,000,000,1
  391.  * Author: Max R. Dürsteler, 9/20/83
  392.  */
  393.  double exp(xi)
  394.  double xi;
  395.  {
  396.     double y, ex, px, nn, ds, in;
  397.  
  398.     if (xi > 88.0) return(1.7014117331926443e38);
  399.     if (xi < -88.0) return(0.0);
  400.     ex = 1.0;
  401.     while( xi > 1.0) {
  402.         ex *= 2.718281828459; /* const. e */
  403.         xi--;
  404.     }
  405.     while( xi < -1.0) {
  406.         ex *= .367879441171;  /* 1/e */
  407.         xi++;
  408.     }
  409. /* Slow, but more precise Taylor expansion series */
  410.     y = ds = 1.0; nn = 0.0;
  411.     while ((ds < 0.0 ? -ds : ds) > 0.000000000001) {
  412.         px = xi/++nn;          /* above precision required */
  413.         ds *= px;
  414.         y += ds;
  415.     }
  416.        y *= ex;
  417.  
  418. /*  Chebyshev polynomials: fast, but less precise then expected!
  419.  *    xi = -xi;
  420.  *    y = (((((.0000006906  * xi
  421.  *        +.0000054302) * xi
  422.  *        +.0001715620) * xi
  423.  *        +.0025913712) * xi
  424.  *        +.0312575832) * xi
  425.  *        +.2499986842) * xi
  426.  *        +1.0;
  427.  *     y = ex / (y * y * y * y);
  428.  */
  429.     return(y);
  430. }
  431.  
  432.  
  433. /* LOG10.C
  434.  * Approximation for logarithm of basis 10
  435.  * Range 0 < x < 1e+38
  436.  * Precision: +/- 0.000,000,1
  437.  * Header: math.h
  438.  * Author: Max R. Dürsteler 9/15/83
  439.  */
  440.  
  441. /* Method of Chebyshev polynomials */
  442. /* C. Hastings, jr. 1955 */
  443.  
  444. double log10(x)
  445. double x;
  446. {
  447.     double xi, y, q, q2;
  448.     int ex;
  449.  
  450.     if (x <= 0.0) return(0.); /* Error!! */
  451.     ex = 0.0; xi = x;
  452.     while (xi < 1.0 ) {
  453.         xi *= 10.0;
  454.         ex--;
  455.     }
  456.     while (xi > 10.0) {
  457.         xi *= 0.1;
  458.         ex++;
  459.     }
  460.     q = (xi - 3.16227766) / (xi + 3.16227766); q2 = q * q;
  461.     y = ((((.191337714 * q2
  462.           + .094376476) * q2
  463.           + .177522071) * q2
  464.           + .289335524) * q2
  465.           + .868591718) * q + .5;
  466.     y += (double) ex;
  467.     return(y);
  468. }
  469.  
  470. /* PW10.C
  471.  * Calculates 10 power x
  472.  * Range: 0 <= x <= 1
  473.  * Precision: +/- 0.000,000,005
  474.  * Header: math.lib
  475.  * Author: Max R. Dürsteler 9/15/83
  476.  */
  477.  
  478. double pw10(x)
  479. double x;
  480. {
  481.     double xi,ex,y;
  482.  
  483.     if (x > 38.0) return(1.7014117331926443e38);
  484.     if (x < -38.0) return(0.0);
  485.     xi = x; ex = 1.0;
  486.     while (xi > 1.0) {
  487.         ex *= 10.0;
  488.         xi--;
  489.     }
  490.     while (xi < 0.0) {
  491.         ex *= 0.1;
  492.         xi++;
  493.     }
  494.     y = ((((((.00093264267    * xi
  495.         + .00255491796) * xi
  496.         + .01742111988) * xi
  497.         + .07295173666) * xi
  498.         + .25439357484) * xi
  499.         + .66273088429) * xi
  500.         +1.15129277603) * xi + 1.0;
  501.     y *= y;
  502.     return(y*ex);
  503. }
  504.  
  505.